Sformułowania silne oraz słabe dla naszego problemu transportu ciepła na obszarze w kształcie litery L rozwiązywanego klasyczną metodą elementów skończonych jest takie same jak w przypadku stosowania metody izogeometrycznej.
\( a(u,v)=l(v) \\ a(u,v) = \int_{\Omega} \left(\frac{\partial u }{\partial x }\frac{\partial v }{\partial x } + \frac{\partial u }{\partial y } \frac{ \partial v}{\partial y } \right)v(x,y) dxdy \\ l(v) = \int_{ \partial \Omega_N } g(x,y) v(x,y) dS +\int_{\Omega } f(x,y) v(x,y) dxdy \)

W naszym przykładowym problemie przepływu ciepła na obszarze w kształcie litery L przyjmujemy \( f(x,y)=0 \) (czyli brak źródeł ciepła).
Na "zewnętrznym" fragmencie brzegu, który nazywamy brzegiem Neumanna \( \Gamma_N = \{(x,y):x=-1 \cup x=1 \cup y=1 \cup y=-1 \} \)
zakładamy znaną prędkość zmiany temperatury.
\( \frac{\partial u(x,y)}{\partial n }=g(x,y) \textrm{ dla } (x,y) in \Gamma_N \)

W celu zdefiniowania warunku brzegowego Neumanna, wprowadzamy biegunowy układ współrzędnych ze środkiem w punkcie (0,0) i z osią odciętą leżącą na osi OX kartezjańskiego układu współrzędnych.

W naszym przykładzie przepływu ciepła na obszarze w kształcie litery L, przyjmujemy przykładową funkcję prędkości przepływu ciepła zdefiniowaną w biegunowym układzie współrzędnych,
\( R \times (0,2\pi) \ni (r,\theta)\rightarrow g(r,\theta)=r^{2/3 }sin^{2/3 }(\theta+\pi/2) \)
gdzie \( r=(x^2+y^2)^{0.5 } \) oznacza odległość od początku układu współrzędnych, natomiast \( \theta = atan2(y,x) \) jest kątem pomiędzy punktem \( (x,y) \) a osią \( x \). Tak więc
\( g(r,\theta)=r^{2/3 }sin^{2/3 }(\theta+\pi/2)=(x^2+y^2)^{1/3 }sin^{2/3 }(atan2(y,x)+\pi/2)=g(x,y) \).

Na "wewnętrznym" fragmencie brzegu, który nazywamy brzegiem Dirichleta \( \Gamma_D = \partial \Omega \ \Gamma_N \) przyjmujemy zerową wartość temperatury
\( u(x,y)=0 \textrm{ dla } (x,y) in \Gamma_D \)

Przykładowa siatka zbudowana z elementów trójkątnych oraz rozpięte na niej wielomiany Lagrange pierwszego stopnia przedstawiona jest na rysunkach Rys. 1-Rys. 2.

Załóżmy, że stosujemy do aproksymacji rozwiązania wielomiany Lagrange'a pierwszego stopnia. W przypadku stosowania liniowych funkcji bazowych są one związane z wierzchołkami elementów trójkątnych (w sensie, że ich maksima znajdują się nad wierzchołkami, a określone są na sąsiadujących elementach trójkątych).
Gdyby zastosować natomiast wielomiany Lagrange'a wyższego stopnia, byłyby one również związane z krawędziami i wnętrzami elementów. Szczegóły podano w rozdziale opisującym sposoby definiowania funkcji bazowych.

W klasycznej metodzie elementów skończonych nasza siatka oraz wynikająca z niej macierz układu równań liniowych nie posiadają regularnego kształtu, tak jak ma to miejsce w przypadku grupy elementów w izogeometrycznej metodzie elementów skończonych. Należy więc wygenerować globalną numerację funkcji bazowych, związanych z wierzchołkami elementów trójkątnych. Każdy element trójkątny posiada z kolei swoją lokalną numeracje wierzchołków. Konieczne jest opracowanie mapy z numeracji globalnej na numerację lokalną każdego elementu trójkątnego. Jest ona narysowana na Rys. 3.

Następnie nad wzorcowym elemetem trójkątnym \( \hat{E}:(0,1)-(0,0)-(1,0) \) wprowadzamy trzy lokalne funkcje bazowe
\( \hat{\psi}^1(\xi,\eta)=(1-\xi)\eta \\ \hat{\psi}^2(\xi,\eta)=(1-\xi)(1-\eta) \\ \hat{\psi}^3(\xi,\eta)=\xi(1-\eta) \)

Funkcje te mapowane są na funkcje bazowe rozpięte na poszczególnych elementach. W tym celu definiujemy mapy odwzorowujące element wzorcowy na poszczególne elementy
\( \hat{E} \ni (\xi,\eta) \rightarrow x_{E1}(\xi,\eta)=(\xi-1,\eta)=(x,y)\in E_1 \\ \hat{E} \ni (\xi,\eta) \rightarrow x_{E2}(\xi,\eta)=(-\xi,-\eta+1)=(x,y)\in E_2 \\ \hat{E} \ni (\xi,\eta) \rightarrow x_{E3}(\xi,\eta)=(\xi,\eta)=(x,y)\ni E_3 \\ \hat{E} \in (\xi,\eta) \rightarrow x_{E4}(\xi,\eta)=(\xi+1,\eta+1)=(x,y)\in E_4 \\ \hat{E} \ni (\xi,\eta) \rightarrow x_{E5}(\xi,\eta)=(\xi,\eta-1)=(x,y)\in E_5 \\ \hat{E} \ni (\xi,\eta) \rightarrow x_{E6}(\xi,\eta)=(1-\xi,-\eta)=(x,y)\in E_6 \)

Następnie obliczamy odwzorowania odwrotne
\( E_1 \ni (x,y) \rightarrow x_{E1}^{-1}(x,y)=(x+1,y)=(\xi,\eta)\in \hat{E} \\ E_2 \ni (x,y) \rightarrow x_{E2}^{-1}(x,y)=(-x,-y+1)=(\xi,\eta)\in \hat{E} \\ E_3 \ni (x,y) \rightarrow x_{E3}^{-1}(x,y)=(x,y)=(\xi,\eta)\in \hat{E} \\ E_4 \ni (x,y) \rightarrow x_{E4}^{-1}(x,y)=(x-1,y-1)=(\xi,\eta) \in \hat{E} \\ E_5 \ni (x,y) \rightarrow x_{E5}^{-1}(x,y)=(x,y+1)=(\xi,\eta) \in \hat{E} \\ E_6 \ni (x,y) \rightarrow x_{E6}^{-1}(x,y)=(-x+1,-y)=(\xi,\eta)\in \hat{E} \)

Mając wzory na odwzorowania mapujące pomiędzy elementem wzorcowym a poszczególnymi elementami, możemy zdefiniować funkcje bazowe na elemencie pierwszym
\( \psi^1(x,y)=\hat{\psi}^1\left(x_{E1}^{-1}(x,y)\right)=\hat{\psi}^1\left(x+1,y \right)=(1-(x+1))y=-xy \\ \psi^2(x,y)=\hat{\psi}^2\left(x_{E1}^{-1}(x,y)\right)=\hat{\psi}^2\left(x+1,y \right)=(1-(x+1))(1-y)=x(y-1) \\ \\ \psi^3(x,y)=\hat{\psi}^3\left(x_{E1}^{-1}(x,y)\right)=\hat{\psi}^3\left(x+1,y \right)=(x+1)(1-y) \)
(zauważmy, że pierwsza lokalna funkcja przyjmuje wartość 1 w punkcie (-1,1), druga lokalna funkcja przyjmuje wartość 1 w punkcie (-1,0), oraz trzecia lokalna funkcja przyjmuje wartość 1 w punkcie (0,0)),

na elemencie drugim
\( \psi^1(x,y)=\hat{\psi}^1\left(x_{E2}^{-1}(x,y)\right)=\hat{\psi}^1\left(-x,-y+1 \right)=(1+x)(-y+1)=(1+x)(1-y) \\ \psi^2(x,y)=\hat{\psi}^2\left(x_{E2}^{-1}(x,y)\right)=\hat{\psi}^2\left(-x,-y+1 \right)=(1+x)(1-(-y+1))=(1+x)y \\ \\ \psi^3(x,y)=\hat{\psi}^3\left(x_{E2}^{-1}(x,y)\right)=\hat{\psi}^3\left(-x,-y+1 \right)=-x(1-(-y+1))=-xy \)
(zauważmy, że pierwsza lokalna funkcja przyjmuje wartość 1 w punkcie (0,0), druga lokalna funkcja przyjmuje wartość 1 w punkcie (0,1), oraz trzecia lokalna funkcja przyjmuje wartość 1 w punkcie (-1,1)),

na elemencie trzecim
\( \psi^1(x,y)=\hat{\psi}^1\left(x_{E3}^{-1}(x,y)\right)=\hat{\psi}^1\left(x,y \right)=(1-x)y \\ \psi^2(x,y)=\hat{\psi}^2\left(x_{E3}^{-1}(x,y)\right)=\hat{\psi}^2\left(x,y \right)=(1-x)(1-y) \\ \\ \psi^3(x,y)=\hat{\psi}^3\left(x_{E3}^{-1}(x,y)\right)=\hat{\psi}^3\left(x,y \right)=x(1-y) \)
(zauważmy, że pierwsza lokalna funkcja przyjmuje wartość 1 w punkcie (0,1), druga lokalna funkcja przyjmuje wartość 1 w punkcie (0,0), oraz trzecia lokalna funkcja przyjmuje wartość 1 w punkcie (1,0)),

na elemencie czwartym
\( \psi^1(x,y)=\hat{\psi}^1\left(x_{E4}^{-1}(x,y)\right)=\hat{\psi}^1\left(x-1,y-1 \right)=(1-(x-1))(y-1)=x(1-y) \\ \psi^2(x,y)=\hat{\psi}^2\left(x_{E4}^{-1}(x,y)\right)=\hat{\psi}^2\left(x-1,y-1 \right)=(1-(x-1))(1-(y-1))=xy \\ \\ \psi^3(x,y)=\hat{\psi}^3\left(x_{E4}^{-1}(x,y)\right)=\hat{\psi}^3\left(x-1,y-1 \right)=(x-1)(1-(y-1))=(1-x)y \)
(zauważmy, że pierwsza lokalna funkcja przyjmuje wartość 1 w punkcie (1,0), druga lokalna funkcja przyjmuje wartość 1 w punkcie (1,1), oraz trzecia lokalna funkcja przyjmuje wartość 1 w punkcie (0,1)),

na elemencie piątym
\( \psi^1(x,y)=\hat{\psi}^1\left(x_{E5}^{-1}(x,y)\right)=\hat{\psi}^1\left(x,y+1 \right)=(1-x)(y+1) \\ \psi^2(x,y)=\hat{\psi}^2\left(x_{E5}^{-1}(x,y)\right)=\hat{\psi}^2\left(x,y+1 \right)=(1-x)(1-(y+1))=(x-1)y \\ \\ \psi^3(x,y)=\hat{\psi}^3\left(x_{E5}^{-1}(x,y)\right)=\hat{\psi}^3\left(x,y+1)\right)=x(1-(y+1))=-xy \)
(zauważmy, że pierwsza lokalna funkcja przyjmuje wartość 1 w punkcie (0,0), druga lokalna funkcja przyjmuje wartość 1 w punkcie (0,-1), oraz trzecia lokalna funkcja przyjmuje wartość 1 w punkcie (1,-1)),

oraz na elemencie szóstym
\( \psi^1(x,y)=\hat{\psi}^1\left(x_{E6}^{-1}(x,y)\right)=\hat{\psi}^1\left(-x+1,-y \right)=(1-(-x+1))(-y)=-xy \\ \psi^2(x,y)=\hat{\psi}^2\left(x_{E6}^{-1}(x,y)\right)=\hat{\psi}^2\left(-x+1,-y \right)=(1-(-x+1))(1-(-y))=x(1-y) \\ \\ \psi^3(x,y)=\hat{\psi}^3\left(x_{E6}^{-1}(x,y)\right)=\hat{\psi}^3\left(-x+1,-y)\right)=(1-x)(1+y) \)
(zauważmy, że pierwsza lokalna funkcja przyjmuje wartość 1 w punkcie (1,-1), druga lokalna funkcja przyjmuje wartość 1 w punkcie (1,0), oraz trzecia lokalna funkcja przyjmuje wartość 1 w punkcie (0,0)).

Funkcje na lokalnych elementach agregowane są w globalne funkcje bazowe \( e^1,e^2,e^3,e^6,e^8 \) zgodnie z mapowaniem przedstawionym na Rys. 3.
Innymi słowy
\( e_1=\psi^3_{E1}+\psi^2_{E2} \\ e_2=\psi^2_{E2}+\psi^3_{E3}+\psi^2_{E4} \\ e_3=\psi^1_{E4}\\ e_6 =\psi^2_{E3}+\psi^3_{E4}+\psi^1_{E6} \\ e_8=\psi^2_{E5}+\psi^3_{E6} \)

Nasze rozwiązanie (rozkład temperatury) przybliżamy za pomocą kombinacji zagregowanych wielomianów Lagrange'a
\( u(x,y)=\sum_{i=1,...,8 } u_i e_i(x,y) \).
Ze względu na stosowanie zerowego warunku Dirichleta na brzegu "wewnętrznym" możemy od razu przyjąć, że współczynniki wielomianów Lagrange'a związane z brzegiem wewnętrznym są równe zero, czyli
\( u_4=u_5=u_7=0 \), dlatego więc nawet nie wyprowadzaliśmy wzorów na funkcje bazowe \( e_4(x,y),e_5(x,y),e_7(x,y) \), ponieważ są one zerowane przez zerowy warunek Dirichleta (wymuszany przez zerowanie współczynników).

Generowanie nieregularnej globalnej numeracji funkcji bazowych oraz nieregularnego mapowania lokalnych funkcji elementowych na funkcje globalne jest typową cechą stosowania metody elementów skończonych na siatkach nieregularnych na przykład z wykorzystaniem wielomianów Lagrange'a. W ogólnym przypadku, musimy również wygenerować globalną numerację i lokalne numerację wierzchołków elementów trójkątnych, oraz wygenerować wzory na funkcje bazowe na wszystkich elementach trójkątnych.